suppressPackageStartupMessages({
library(readr)
library(tidyverse)
library(reshape)
library(popbio)
})
Our data is composed by length data, which we transformed into expected age using the transformed Von Bertalanffy function \[ A = \frac {1}{-K} \times log(1-(\frac {L}{L_{inf}}))+t_0 \], where: A = age, K = growth rate, L = length, \(L_{inf}\) = asymptotic length. Therefore, our model use a age structure.
Mortality: \(Z = 1.69\) (Garbin & Castello, 2014)
Fecundity: \(a = -1.33354\), \(b = 3.238\) (from FISHBASE)
von Bertanalnffy growth parameters: \(L_inf = 102.0\), \(K = 0.55\), \(t_0 = -0.02\) (Uchiyama & Struhsaker, 1981)
DiagrammeR::grViz("
digraph boxes_and_circles{
# Define nodes
node [shape = circle
penwidth = 2,
fontsize = 24]
egg; a1; a2; a3; a4; a5; a6; a7; a8; a9; a10; a11; a12; a13
# Add edge statements
#Advaance classes
egg -> a1[label = s0]
a1 -> a2[label = s1]
a2 -> a3[label = s2]
a3 -> a4[label = s3]
a4 -> a5[label = s4]
a5 -> a6[label = s5]
a6 -> a7[label = s6]
a7 -> a8[label = s7]
a8 -> a9[label = s8]
a9 -> a10[label = s9]
a10 -> a11[label = s10]
a11 -> a12[label = s11]
a12 -> a13[label = s12]
# Fecundities
a1 -> egg[label = f1]
a2 -> egg[label = f2]
a3 -> egg[label = f3]
a4 -> egg[label = f4]
a5 -> egg[label = f5]
a6 -> egg[label = f6]
a7 -> egg[label = f7]
a8 -> egg[label = f8]
a9 -> egg[label = f9]
a10 -> egg[label = f10]
a11 -> egg[label = f11]
a12 -> egg[label = f12]
a13 -> egg[label = f13]
}
", height = 1000)
Figure 1 - Life-cycle diagram for Skipjack tuna (Katsuwonus pelamis). S indicates survival, and f indicates fecundity.
The model we have defined axplicitly models eggs as a cohort (age = 0). We will conssider deleting this in the future, and explicitly calculate recruits.
#Build a matrix with text
A <- matrix(0, 14, 14) #initial empty matrix with all 0
# By using the $$, we can embed latex equations that will then be rendered to loog better
# Populate matrix with mortality
for (i in 2:14){
A[i,i-1] <- paste0("$s_{", i-1, ",", i-2, "}$")
}
ages <- seq(0:13)-1
A[1,] <- paste0("$f_{", seq(1:14)-1,"}$")
knitr::kable(A, col.names = paste0("$a_{",ages,"}$"), row.names = F, caption = "Table 1- Symbolical population matrix A. Subindices represent their rows and columns, respectively. The inferior diagonal represents survivals, while the first row represents facundities.")
| \(a_{0}\) | \(a_{1}\) | \(a_{2}\) | \(a_{3}\) | \(a_{4}\) | \(a_{5}\) | \(a_{6}\) | \(a_{7}\) | \(a_{8}\) | \(a_{9}\) | \(a_{10}\) | \(a_{11}\) | \(a_{12}\) | \(a_{13}\) |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| \(f_{0}\) | \(f_{1}\) | \(f_{2}\) | \(f_{3}\) | \(f_{4}\) | \(f_{5}\) | \(f_{6}\) | \(f_{7}\) | \(f_{8}\) | \(f_{9}\) | \(f_{10}\) | \(f_{11}\) | \(f_{12}\) | \(f_{13}\) |
| \(s_{1,0}\) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 0 | \(s_{2,1}\) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 0 | 0 | \(s_{3,2}\) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 0 | 0 | 0 | \(s_{4,3}\) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 0 | 0 | 0 | 0 | \(s_{5,4}\) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 0 | 0 | 0 | 0 | 0 | \(s_{6,5}\) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 0 | 0 | 0 | 0 | 0 | 0 | \(s_{7,6}\) | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 0 | 0 | 0 | 0 | 0 | 0 | 0 | \(s_{8,7}\) | 0 | 0 | 0 | 0 | 0 | 0 |
| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | \(s_{9,8}\) | 0 | 0 | 0 | 0 | 0 |
| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | \(s_{10,9}\) | 0 | 0 | 0 | 0 |
| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | \(s_{11,10}\) | 0 | 0 | 0 |
| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | \(s_{12,11}\) | 0 | 0 |
| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | \(s_{13,12}\) | 0 |
# Define Von Bert parameters
# Garbin & Castello, 2014
# L_inf <- mean(c(80, 80, 87.12, 94, 97.9, 97.3, 112.34, 89.38, 92.5))
# K <- mean(c(0.32, 0.6, 0.22, 0.38, 0.14, 0.25, 0.14, 0.38, 0.16))
# Or from Us and Stiurskjgfas, 1981
L_inf <- 102
K <- 0.55
t_0 <- -0.02
# Just to be clear, we use L-inf and to from Ushisomething, and K from the review
# Define fecundity parameters
# From fishbase http://www.fishbase.se/Reproduction/FecundityList.php?ID=107&GenusName=Katsuwonus&SpeciesName=pelamis&fc=416&StockCode=121
fec_a <- -1.33354
fec_b <- 3.238
# Define mortality
m <- 0.63
z <- 1.39
# Convert length to age using von bertalanffy model, solving for t
length2age <- function(length, l_inf, K, t_o){
age <- (1/-K)*(log(1-(length/L_inf))) + t_o
return(age)
}
# Convert age to length using von bertalanffy model
age2length <- function(age, l_inf, K, t_o){
length <- l_inf*(1-exp(-K*(age-t_o)))
return(length)
}
#Convert length to fecundity (number of eggs)
fecundity <- function(length, a, b){
f <- 10^(a+(b*log10(length*10)))
return(f)
}
A <- matrix(0, 14, 14) #initial empty matrix with all 0
# Populate matrix with mortality
for (i in 2:14){
A[i,i-1] <- exp(-1-z)
}
# Populate matrix with fecundity
ages <- seq(0:13)-1
lengths <- age2length(ages, L_inf, K, t_0)
A[1,] <- fecundity(lengths, fec_a, fec_b)
A[is.na(A)] <- 0
A[2,1] <- 0.0000001
A[1,1] <- 0
colnames(A) <- ages
rownames(A) <- ages
knitr::kable(A, col.names = paste0("$a_{",ages,"}$"), row.names = F, caption = "Table 2- Population matrix A. The inferior diagonal represents survivals, while the first row represents facundities.")
| \(a_{0}\) | \(a_{1}\) | \(a_{2}\) | \(a_{3}\) | \(a_{4}\) | \(a_{5}\) | \(a_{6}\) | \(a_{7}\) | \(a_{8}\) | \(a_{9}\) | \(a_{10}\) | \(a_{11}\) | \(a_{12}\) | \(a_{13}\) |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0e+00 | 1.657256e+07 | 7.026737e+07 | 1.294406e+08 | 1.758241e+08 | 2.072322e+08 | 2.27012e+08 | 2.38998e+08 | 2.461086e+08 | 2.502768e+08 | 2.527038e+08 | 2.541114e+08 | 2.549259e+08 | 255396708 |
| 1e-07 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.00000e+00 | 0.00000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 9.162970e-02 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.00000e+00 | 0.00000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 9.162970e-02 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.00000e+00 | 0.00000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 9.162970e-02 | 0.000000e+00 | 0.000000e+00 | 0.00000e+00 | 0.00000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 9.162970e-02 | 0.000000e+00 | 0.00000e+00 | 0.00000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 9.162970e-02 | 0.00000e+00 | 0.00000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 9.16297e-02 | 0.00000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.00000e+00 | 9.16297e-02 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.00000e+00 | 0.00000e+00 | 9.162970e-02 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.00000e+00 | 0.00000e+00 | 0.000000e+00 | 9.162970e-02 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.00000e+00 | 0.00000e+00 | 0.000000e+00 | 0.000000e+00 | 9.162970e-02 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.00000e+00 | 0.00000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 9.162970e-02 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.00000e+00 | 0.00000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 9.162970e-02 | 0 |
load("size_dist_1980.RData") #This is the n-vector for 1980, loads it to a vector called n
n <- n %>% {
.$N}%>%
c(rep(0, times = 8))
n_0 <- sum(A[1, 2:14]*n, na.rm = T)
n <- c(n_0, n)
n
## [1] 2.562588e+19 8.674537e+06 3.136387e+11 2.711014e+10 4.391505e+08
## [6] 3.880114e+06 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [11] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
project <- popbio::pop.projection(A, n, 15)
pop <- project$stage.vectors %>%
as.data.frame() %>%
mutate(Age = as.factor(seq(0:13)-1)) %>%
gather(Year, N, -Age) %>%
mutate(Year = as.numeric(as.character(Year))) %>%
select(Year, Age, N)
ggplot(pop, aes(x = Year, y = log(N), color = Age)) +
geom_line() +
theme_bw()
Figure 2 - Population size through time, represented by ages.
\(\lambda =\) 1.4729341
plot(project$pop.changes, type = "b", xlab = "Iterations", ylab = "Lambda")
Figure 3 - Convergence of \(\lambda\) through time.
pop %>%
filter(Year == 14) %>%
ggplot(aes(x = Age, y = log(N))) +
geom_bar(stat = "identity") +
theme_bw()
Figure 4 - Stabe age structure, reached after 15 iterations in the model.
sensitivity(A) %>%
image2(box.offset=.1)
Figure 5 - Sensitivity matrix.
elasticity(A) %>%
image2(box.offset = 0.1)
Figure 6 - Elasticity matrix.
The \(\lambda = 1.47\) shows that the population is projected to grow at rate of 47%. The result of the sensitivity analysis suggest that is most effective to intervene on the early ages. Since the species mature early (between ages 1 and 2), there is no necessity of conserving the old/big individuals that have higher fecundity, because the younger ones (in quantity) can provide enough eggs to maintain the population.
The model assumes the mortality is constant through ages. Better estimates of mortality/survival for different ages would show the real population structure.
The model assumes that the population is density-independent. However, we believe that fecundity, recruitment and survival of this pelagic species is density dependent. Perhaps a first approach wil be to use a Beverton-Holt approach to calculate recruitment, although this assumes density dependence in all the model.
Currently, all organisms are harvested in the model (excluding eggs). An interesting approach would be to set a size limit that reduces total effort on small sizes to natural effort by excluding fishing effort. This could allow us to evaluate the effect of a minimum catch size, management intervention aligned with the observed elasticities.
Environmental variation is another aspect to consider. It is known that recruitment of fish can be affected by extreme environmental conditions, like El Nino. Further analysis could incorporate the effects of climate variability on recruitment.